1.1 Principle and prerequisites

The BRD normalization aims to rescale read counts or genomic coverage values from different ChIP-seq experiments, such that background levels - which should correspond mainly to ChIP-seq reads from non-specific DNA fragments - become equalized between the experimental conditions.

To achieve this goal the BRD method is designed upon the following assumptions and prerequisites.

1. For all considered experimental conditions, some portions of the genome are invariably devoid of the immuno-precipitated chromatin mark.

   We use the term background candidates to designate such invariable genomic regions in a given experimental setup.

2. Background candidates can be distinguished as a specific density mode in appropriately transformed ChIP-seq read count distributions.

A first prerequisite when applying the BRD method is to define the portions of the genome that are biologically relevant for the search of background candidates.

Although a completely naive search of background candidates among the whole genome should be possible in theory, without further assumptions this search may not lead to the selection of valid candidates. For instance, some chromatin marks can cover most of the genome at low but distinguishable levels, as it can be observed with H3K27me3 ChIP-seq profiles in populations of Embryonic Stem Cells (ESCs) or Neural Stem Cells (NSCs). In such situation the second assumption (b. above) of the BRD method does not hold since the density mode of background candidates can be masked by a close and more prevalent mode corresponding to minimal but extensive ChIP enrichments.

For commonly studied marks that are linked with the regulation of gene expression, a reasonable option consist in searching background candidates over gene units or transcription start sites (TSS), assuming that the mark of interest would be always absent at a subpopulation of genes, for instance the constitutively expressed ones or silent ones in the considered experimental conditions.

The second prerequisite is that the ChIP-seq dataset must include background profiles (e.g. input DNA, ChIP-seq performed with IgG or in KO conditions), preferably generated using the same sequencing setup as for the ChIP profiles. For technical reasons, 2 background profiles is the minimum for the current version of Tightrope.

1.2 Technical overview

The figure below illustrates a ChIP-seq normalization workflow with the current version of Tightrope, giving an overview of algorithms implemented in its core functions (A), as well as examples of resulting control graphs (B,C) which provide guidance for the choice of parameters required by these functions.

As shown in this workflow, the BRD method consists mainly in the transformation of read counts over chosen genomic regions (A) into a bivariate density distribution (B), a technique that we name Count Density after Dithering and Dimensionality Reduction (CDaDaDR) which combines a principal component analysis and a non-parametric density estimation, followed by the identification of the local density mode which corresponds to the best background candidates (C) thanks to a clustering technique based on density gradients.

1.3 User interventions

Overall, using the BRD method involves two key interventions. As explained in the previous section, the first one consists in choosing genomic regions that are relevant for the search of background candidates and generating the corresponding read count matrix.

The second intervention consists in tuning the most influencial parameters of the BRD function, which are the number of distinct density modes (i.e. clusters) and the density threshold values used to define these clusters. Practically, the default threshold values should perform reasonably well for common use cases and the remaining decisions after generating read count matrixes will often be limited to a choice between 1 or 2 clusters. This choice should be made after observing the bivariate density distribution generated by the PlotBRD function, opting for a value of 1 or 2 clusters for a globally unimodal or bimodal density respectively.

2.1 R and RStudio environments

If R and RStudio are not already installed.

• Manually download and install R.

• Manually download and install RStudio.

2.2 R packages

Once R and RStudio are installed, open RStudio and run the R code below which should install all required packages automatically. While R will perfom these tasks, you may be prompted for a confirmation of package installations or updates.

# Detect already installed packages
pkg <- installed.packages()[, "Package"]

# CRAN packages
lst <- setdiff(c("devtools"), pkg)
if(length(lst) > 0) install.packages(lst)

# Bioconductor packages
source("https://bioconductor.org/biocLite.R")
lst <- c(
  "Rsamtools", "GenomicAlignments", "GenomicRanges", "GenomicFeatures",
  "rtracklayer", "biomaRt"
)
lst <- setdiff(lst, pkg)
if(length(lst) > 0) biocLite(lst)

# GitHub packages
library("devtools")
install_github("benja0x40/Barbouille")
install_github("benja0x40/Tightrope")

2.3 Tightrope

Once Tightrope is installed, the first thing to do when starting a work session is to load the package into the active R session.

 # Load Tightrope package (this will produce numerous messages)
library(Tightrope)

At this stage, it is possible to obtain the list of R functions available in Tightrope and to browse the builtin documentation by calling the help() function in the R console.

help(package = "Tightrope")

And the present document can be shown by using the vignette() function.

vignette("BRD", package = "Tightrope")

3.1 Auto saving and loading of R data

3.2 Genomic annotations and relevant regions

See the dedicated vignette.

vignette("Annotations", package = "Tightrope")

# Load mouse genome annotations (mm10)
data("EGA91_mouse")    # Gene features from Ensembl      (environment)
data("CGI_mm10")       # CpG Islands from UCSC           (GRanges)
data("mm10_blacklist") # Blacklisted regions from ENCODE (GRanges)

# Load whole genome tiling object WGT or create it and save it if necessary
MakeObj(WGT, {
  # Make 4kb windows every 1kb over the whole genome (GRanges object)
  WGT <- GenomicTiling(EGA91_mouse$seqinfo, s = 1000, w = 4000)
  
  WGT <- CleanupGRanges(
    WGT, EGA91_mouse$seqinfo, "Mus musculus", mm10_blacklist
  )
  
  # Functional annotation of genomic bins
  grp <- factor(
    rep(1, length(WGT)), levels = 1:5, ordered = T,
    labels = c("intergenic", "genic", "TSS", "CGI", "blacklist")
  )
  grp[which(countOverlaps(WGT, EGA91_mouse$GNU) > 0)] <- "genic"
  grp[which(countOverlaps(WGT, EGA91_mouse$TSS) > 0)] <- "TSS"
  grp[which(countOverlaps(WGT, CGI_mm10) > 0)] <- "CGI"
  grp[which(countOverlaps(WGT, mm10_blacklist) > 0)] <- "blacklist"
  
  WGT$category <- grp
})

3.3 Read count data

# Define bam file paths and the corresponding names of experimental conditions
bam.files <- c(
  H3WT_DMSO_H3K27me3_R1 = "~/sequencing/sample1.bam",
  K27M_DMSO_H3K27me3_R1 = "~/sequencing/sample2.bam",
  H3WT_DMSO_Input_R1    = "~/sequencing/sample3.bam",
  K27M_DMSO_Input_R1    = "~/sequencing/sample4.bam"
  # etc.
)
conditions <- names(bam.files)

# Define accepted reads (single-end sequencing)
bam.flag <- scanBamFlag(isDuplicate = F, isUnmappedQuery = F)

# Read counts over genomic tiling bins WGT (single-end sequencing)
COUNTS <- ReadCountMatrix(
    bam.files, WGT, paired = F, names = conditions,
    param = ScanBamParam(flag = bam.flag)
)

# Define accepted reads (paired-end sequencing)
bam.flag <- scanBamFlag(
  isDuplicate = F, isUnmappedQuery = F, isProperPair = T
)
  
# Read counts over genomic tiling bins WGT (paired-end sequencing)
COUNTS <- ReadCountMatrix(
    bam.files, WGT, paired = T, names = conditions,
    param = ScanBamParam(flag = bam.flag)
)

# Load the NSC dataset (replicate of experiments from Mohammad et al. 2017)
data("NSC_K27M")
library(data.table)

NSC_K27M$METADATA # Annotation of experimental conditions (data.table object)
NSC_K27M$SPIKEIN  # Number of Drosophila reads (chromatin spike-in)
NSC_K27M$COUNTS   # Precomputed read count matrixes over genomic tiling bins
NSC_K27M$WGT      # Whole genome tiling bins (GRanges object)

COUNTS <- NSC_K27M$COUNTS
WGT    <- NSC_K27M$WGT

3.4 Genomic read count distributions

Here we select the chip and ctrl sample identifiers corresponding to H3K27me3 and Input ChIP-seq profiles in duplicate experiments.

chip <- NSC_K27M$METADATA[antibody == "H3K27me3"]$sample_id
ctrl <- NSC_K27M$METADATA[antibody == "Input"]$sample_id

# log2 transformed counts with dithering
MakeObj(l2c, {
  l2c <- log2(DitherCounts(COUNTS))
  l2c <- RemoveInputBias(l2c, controls = ctrl, as.log2 = T)
})

## [LittleThumb] bypass | l2c = AutoSaved/l2c.rds

# RPM
MakeObj(rpm, {
  rpm <- log2(colSums(COUNTS))
  rpm <- t(t(l2c) - rpm + mean(rpm))
})

## [LittleThumb] bypass | rpm = AutoSaved/rpm.rds

PlotCountDistributions(
  rpm[, chip], pops = WGT$category, sampling = 5E5, main = "Genomic bins"
)
chk <- FiniteValues(rpm[, ctrl])
abline(h = median(rpm[chk, ctrl]), lwd = 1.5, lty = 3)

3.5 Preview of BRD using genes, promoters and CGI bins

# Prepare figure layout
layout(matrix(1:6, 2, 3, byrow = F))

# Run BRD with default parameters and plot the density graph
PlotBRD(
  BRD(cnt = COUNTS[WGT$category == "genic", ], controls = ctrl, sampling = 1E5),
  plots = c("density", "thresholds"), title = "using genes"
)
PlotBRD(
  BRD(cnt = COUNTS[WGT$category == "TSS", ], controls = ctrl, sampling = 1E5),
  plots = c("density", "thresholds"), title = "using promoters"
)
PlotBRD(
  BRD(cnt = COUNTS[WGT$category == "CGI", ], controls = ctrl, sampling = 1E5),
  plots = c("density", "thresholds"), title = "using CpG-islands"
)

When using gene units, we observe that the bivariate density distribution (upper left panel) is bimodal, meaning that BRD should account for 2 density clusters when searching for the background candidates (ncl = 2).

When using promoters, we observe that the bivariate density distribution (upper left panel) is unimodal, meaning that BRD should account for a single density cluster when searching for the background candidates (ncl = 1).

When using CGIs, we observe that the bivariate density distribution (upper left panel) is bimodal, meaning that BRD should account for 2 density clusters when searching for the background candidates (ncl = 2).

3.6 Estimation of scaling factors using genes, promoters and CGI bins

# Run BRD using read counts over gene units (unimodal density)
brd.gnu <- BRD(
  cnt = COUNTS[WGT$category == "genic", ], controls = ctrl, ncl = 1,
  sampling = 3E5, dither = 10
)
# Run BRD using read counts over promoters (bimodal density)
brd.tss <- BRD(
  cnt = COUNTS[WGT$category == "TSS", ], controls = ctrl, ncl = 1,
  sampling = 3E5, dither = 10
)
# Run BRD using read counts over CpG-islands (unimodal density)
brd.cgi <- BRD(
  cnt = COUNTS[WGT$category == "CGI", ], controls = ctrl, ncl = 1,
  sampling = 3E5, dither = 10
)

# Value of BRD scaling factors
ScalingFactors(brd.gnu) # Estimated using gene units
ScalingFactors(brd.tss) # Estimated using promoters
ScalingFactors(brd.cgi) # Estimated using CpG-islands

The scaling factors estimated from either gene units, promoters or CpG-islands are globally similar, with differences below 10%. For normalizations, we choose to use scaling factors estimated with CpG-islands read counts, which present intermediate values compared to the estimations from genes or promoter read counts.

The next page summarizes the BRD control graphs.

# Prepare figure layout
layout(matrix(1:12, 4, 3, byrow = F))

# Plot BRD control graphs using genes, promoters and CpG-islands
PlotBRD(brd.gnu, title = "using genes")       # Left panels
PlotBRD(brd.tss, title = "using promoters")   # Center panels
PlotBRD(brd.cgi, title = "using CpG-islands") # Right panels

3.7 Normalization of read counts

# l2c <- log2(DitherCounts(COUNTS))
# 
# chk <- NormalizeCountMatrix(l2c, brd.cgi, as.log2 = T)
# nrm <- RemoveInputBias(chk, controls = ctrl, as.log2 = T)
# 
# rmb <- RemoveInputBias(l2c, controls = ctrl, as.log2 = T)
# rmb <- NormalizeCountMatrix(rmb, brd.cgi, as.log2 = T)
# 
# layout(matrix(1:6, 2, 3, byrow = T))
# trash <- SideBySideDensity(
#   chk[WGT$category == "CGI", ], mapper = cmf, las = 2, x.labels = F, main = "chk"
# )
# trash <- SideBySideDensity(
#   nrm[WGT$category == "CGI", ], mapper = cmf, las = 2, x.labels = F, main = "nrm"
# )
# trash <- SideBySideDensity(
#   rmb[WGT$category == "CGI", ], mapper = cmf, las = 2, x.labels = F, main = "rmb"
# )

# Normalize genomic read counts using BRD scaling factors from CGI bins
MakeObj(nrm, {
  nrm <- NormalizeCountMatrix(l2c, brd.cgi, as.log2 = T)
})

## [LittleThumb] bypass | nrm = AutoSaved/nrm.rds

# Show genomic read count distributions after BRD normalization
PlotCountDistributions(
  nrm[, chip], pops = WGT$category, sampling = 5E5,
  main = "Genomic bins"
)

The distribution of normalized read counts over CpG-islands looks consistent. However, as opposed to the expected effect of K27M on global H3K27me3 levels, we observe that in these experiments, ChIP-seq enrichments in the K27M condition are globally higher than in the WT condition.

3.8 WORK IN PROGRESS

3.9 Adjustment of parameters

Once the read count matrix aimed for the search of background candidates has been generated, a preliminary application of the BRD method takes the following form.

# Define column identifiers corresponding to Input profiles
ctrl <- c("Input_WT", "Input_K27M")

# Preliminary application of the BRD method based on the read count matrix cnt
brd <- BRD(cnt, controls = ctrl)

After this, calling the PlotBRD function will generate control graphs that are necessary to adjust key parameters of the BRD function, that is the number of clusters ncl and the density thresholds bdt.

PlotBRD(brd) # Generate BRD control graphs (e.g. bivariate density distribution)

For instance, with a bivariate density distribution similar to the one shown in section 1.2, the number of clusters must be set to 2 for an accurate estimation of normalization factors.

# Adjusted application of the BRD method for a bimodal bivariate density  
brd <- BRD(cnt, controls = ctrl, ncl = 2)

In this virtual example, the estimated normalization factor for each experimental condition would be given by the object brd$normfactors.

Applying the BRD normalization to log2 transformed read counts or coverage vectors simply consist in adding these factors. Thus, the raw read counts should be normalized as follows.

# Apply BRD normalization factors to raw read counts
nrm <- t(t(cnt) * 2^brd$normfactors)

More detailed examples of normalization with real ChIP-seq datasets can be found in the next section.

3.10 Normalization

3.11 Verifications

4.1 Examples

4.2 Questions & answers

• Can I apply BRD “blindly” for validated marks (or automatically on many datasets)

Short answer: no you should not unless you know that the extact same protocol has been used and the data quality is consistenly good. Because variation of protocols or experimental setup may break BRD assumptions, one needs to verify that the background candidates and read count distributions make sense (look convincing). BRD remains an empirical approach.

• Can I normalize ChIP-seq for different marks all together.

Although this is technically possible, the BRD method was not designed for this and there are several reasons not to do so. The PCA projection does not allow more than one enrichment dimension. The antibody response may yield to different dynamic ranges of enrichment and backgrounds. The background candidates may be mutually exclusive for the distinct marks. In summary: bad idea.

• If I have replicates, should I normalize them separately or all together?

Yes, as much as you like! BRD does take into account replicate variations thanks to the bivariate PCA projection and the estimations are improved.

• Should I have replicates.

Absolutely! We recommend having two replicates for any study aiming at quatitative differences. The more replicates the better. Same for inputs/controls although beyond a handful is going to be an overkill.

• I don’t have replicate inputs.

• I would like to load and use my own annotation object, how should I do?

• Why BRD results differ when runned several times with same parameters?

• Why estimated scaling factors differ slightly whether ncl is unspecified or manually set to 1?

• I have spike-in data showing strong discrepancy with BRD scaling factors